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ABSTRACT 


Investigations of scattered transient waveforms from conducting bodies have 
shown that it is possible to classify electromagnetic scatterers. The concept is based 
upon the natural resonance modes which are part of the scatterer response to an 
incident excitation. A new approach for describing natural resonance modes using 
recursive systems is introduced. A discrete auto-regressive moving-average (ARMA) 
type model for the case of the space-time wave equation is presented. This model 
results from a finite-difference approximation to the wave-equation. The ARMA 
model has spatially-independent coefficients for the temporal recursive terms. 
Computed results showing aspect- and spatial-independence of natural resonance 
modes, with verification of the ARMA model, are also included. Applications to 
target identification, using the natural resonant frequencies of the target’s echo 
signature, are considered. 
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I. INTRODUCTION 


A. BACKGROUND 

Investigations of transient scattered waveforms from conducting bodies have 
shown that it is possible to classify electromagnetic scatterers. Such research is 
particularly applicable to inverse scattering and radar target identification. 
Conceptual applications in radar target identification have been demonstrated in 
some laboratories using advanced, high resolution radar techniques. One type of 
technique is based upon the natural resonance modes which are part of the scattering 
response to an incident transient excitation. In 1971, Baum [Ref, 1], introduced the 
development of the singularity expansion method (SEM), which presents the response 
of a system as a weighted expansion of complex natural modes. Theoretical studies 
and experiments have shown that these modes are functions of the scatterer geometry 
and composition but are independent of the incident excitation, including aspect and 
polarization. 

The fact that these natural modes are only functions of the target led to the 
idea to use them as a data base for the target identification process. This concept 
was introduced by Moffatt and Mains in 1974 [Ref. 2]. The identification process, 
in its most elementary form, includes a comparison against the modes that have been 
extracted and identified, using advanced signal processing techniques, as applied to 
the target’s time-domain scattering response. Morgan showed that a complete 
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description of the scattered signal, using the conventional SEM approach, is valid only 
for the "late-time" portion of the scattered field [Ref 3]. This late-time scattering 
response is due to the source-free currents that remain after the incident field has 
completed its illumination of the target. 

The unique properties of the late-time response are crucial to the development 
of algorithms to identify the resonance modes of the target. These natural resonance 
modes can be represented by pairs of conjugate poles in the left-half complex s-plane 
of the Laplace transform domain. Targets of different geometry', or composition, 
have different pole representations. Several techniques have been developed to 
extract the dominant poles in the time domain scattering responses of simple targets. 
Morgan [Ref 4], for example, introduced the classification of some kinds of 
electromagnetic scatterers by the annihilation of the natural modes. The advantages 
of this technique over others were achieved, primarily, by using only the late-time 
scattered field. The theory of natural resonance scattering is the basis for this thesis, 
and it is therefore explained in more detail in Chapter II. 

B. PRESENTATION IN A DISCRETE MODEL 

Based on the theory' of the natural resonance scattering one can recognize the 
late-time portion of the scattered signal as the response of a linear time-invariant 
(LTI) system. These kinds of systems can be numerically modeled and described by- 
means of linear constant-coefficient difference equations [Ref. 5]. A discrete model, 
which describes the source-free current distribution, may explicitly present the late- 
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time response of the scatterer to an incident field as a function of its natural modes. 
Moreover, the currents may be described by an auto-regressive moving-average 
(ARMA) type model, having constant coefficients for the recursive terms which 
generate the natural resonance modes. This type of model represents an important 
class of discrete systems which are known as recursive systems since the output 
depends on previous values of the output. Such a discrete model can be employed 
in the process of identifying the resonance modes for electromagnetic and acoustic 
scatterers of various shapes. The time-independent coefficients of the recursive terms 
in the ARMA model are also the denominator polynomial coefficients of the system 
transfer function. Using the well-known concept of zeros and poles which represent 
the system frequency response, the roots of this polynomial are the poles which 
represent the complex natural frequencies of the resonance modes [Ref. 6]. In 
addition, such a discrete model gives a complete description of the natural resonance 
modes for a given scattering object when limited by the sampling frequency imposed 
by the Nyquist theorem. This is the basic idea of this thesis, in which such a model 
is set up for the case of the numerical finite-difference solution of the space-time 
wave equation. It is believed that this is the first time that such an approach has 
been presented. 

C. THESIS GOAL 

The goal of this thesis was to investigate the possibilities of describing the 
natural resonance modes of electromagnetic and acoustic scatterers by discrete 
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ARMA models. These ARMA models should have constant coefficients for the 


recursive terms which determine the natural resonance modes. A model for a one¬ 
dimensional structure, such as the damped string with forcing function, illustrates the 
validity of the proposed approach. In addition, three-dimensional scattering 
structures can be analyzed using a space-time finite-difference method, which is an 
extension to the considered herein. 

In this thesis an attempt have been made to present an overview of some 
approaches for solving the problem of the natural modes, and constructing the 
required model. A brief presentation of the theory of natural resonance modes, as 
developed through the time-domain integral equation, is included in Chapter II. 
Chapter III describes the simplified numerical solution, via the vector potential, as 
chosen to demonstrate the possibility of using an ARMA model for the case of an 
electromagnetic thin-wire. This technique is applied to the analogous problem of a 
finite string. Chapter IV presents the ARMA model with several analytical and 
computed results. Conclusions and a description of some future questions are 
discussed in Chapter V. 
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II. THE THEORY OF NATURAL RESONANCE SCATTERING 


A. INTRODUCTION 

The general topic to be considered in this chapter is the natural resonance 
mode representation for the induced current and the scattered field transient 
response of a perfectly conducting body. The case of a thin-wire is presented in some 
detail. The goal of this chapter is to provide the theoretical background required for 
constructing a discrete ARMA model for the natural resonances in electromagnetic 
scattering. In addition, the features of such a model should agree with this theory. 
A description of time-domain solutions relevant to this work are included in this 
chapter. There are two main steps in the process of setting up the ARMA model. 
First, the time-domain numerical solution is required, and second, based on the 
numerical solution, the discrete model must be constructed. 

There are two independent techniques available for solving the problem of 
transient scattering. The first involves the computation of the frequency-domain 
response of the structure, followed by inverse Fourier transformation to yield the 
time-domain response. An alternative approach involves the direct formulation of 
either partial differential or integral equations in the time domain. One way of 
describing the scattered signal in terms of natural modes is by using the formalism 
known as the Singularity Expansion Method (SEM) [Ref. 1]. Mittra and Van 
Blaricum showed that the SEM pole singularities of a structure can be estimated 






directly from its time-domain response [Ref. 7]. In this chapter time-domain integral 
equations are used general to represent solution of the transient electromagnetic 
problem. Also considered will be the conceptual basis of mode representations for 
fields and currents in transient scattering. 


B. TIME DOMAIN INTEGRAL EQUATIONS 

Consider the general three-dimensional transient electromagnetic scattering 
problem as depicted by Fig. 1. 



INCIDENT FIELD 


PERFECTLY CONDUCTING 
OBJECT 


Figure 1. Transient Electromagnetic Scattering [After Ref. 4] 

The perfectly conducting object is illuminated by an incident impulsive type plane 
wave in free space. The incident electric and magnetic fields are E lnc and II lnc , 
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respectively. The integral equation approach represents the induced current on the 
surface of the object in terms of the incident field. There are two fundamental types 
of time-domain integral equations: the Electric Field Integral Equation (EFIE), and 
the Magnetic Field Integral Equation (MFIE). Derivations of both EFIE and MFIE 
are described by Mittra in detail [Ref. 8]. 

Integral equation derivations begin with the time-dependent forms of Maxwell’s 
equations in free space. The continuity equation is used to relate the current density, 
J s , and the charge density a. The scalar potential, <f>, and the vector potential, A, are 
defined in terms of the electric and magnetic fields, respectively. The potentials are 
related to each other via the Lorentz condition[Ref. 8]. The wave equation can be 
derived independently for the scalar potential and for the vector potential using 
Maxwell’s equations. The sources of the nonhomogeneous wave equations for the 
vector and the scalar potentials are the current density, J s , and the charge density, a, 
respectively. The solution is constructed using the Green’s function, which is the 
impulse response in time and space. The integral representation of the time- 
dependent electric and magnetic fields are obtained by applying the solutions of the 
potentials for the general current and charge distribution. 

The expressions of the EFIE and MFIE, for the scattering problem are finally 
derived by applying the appropriate boundary conditions of the tangential electric and 
magnetic fields on the surface of the perfect conductor. In the case of the EFIE the 
sum of the tangential scattered electric field and the incident electric field is zero. 
The boundary condition in the case of the MFIE is that the total tangential magnetic 
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field is equal to the surface current density. The following equations are obtained for 
the EFIE (1) and MFIE (2), [Ref. 8] 



VoJ_ 
R dx 


J M (r',x)- 


a(r',x) R 
e 0 R 3 


13 w N fi 
-—o(r,x )-—- 

cR 2 


e dx 

O 


ds' 


a) 


J s (r,t)=2nxH iHC + 


_L_ A 

27t 


1 fe 


dx 




( 2 ) 


where S is the surface of the conductor, the unit vector which is outward normal to 
S is ft , the current density is J s , the charge density is a, the permeability and 
permittivity of free space are n 0 and e 0 , respectively. Further, the position vector 
is r, r' eS indicates points on the PEC surface, c= l/(s 0 ji 0 ) 1/2 is the velocity of light, 
r=t-R/c is the retarded time, and R=[Rl=ir-H. 

Both the EFIE and MFIE are Principal Value (P.V.) integrals because of kernel 
singularities for r-*r'. The special forms of the time-domain integral equations (1) 
and (2) play fundamental roles in their solution construction. The main difference 
between these equations and their respective frequency domain equations is in the 
solution construction. The equations in frequency domain are handled numerically 
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by matrix inversion, while initial-value techniques are applied to time domain integral 
equations. 

The solution of the general integral equation for the scattering problem uses a 
time-stepping technique. Special cases such as one- and two-dimensional scatterers, 
i.e, cylindrical and wire-type scatterers, have special forms of integral equations, 
hence construction of their solutions differ, considering numerical aspects and 
accuracy. For three-dimensional scatterers the MFIE is most convenient, while the 
EFIE is used in the case of thin wires and thin surfaces. This is due to the fact that 
for solid surface structures the kernel in the MFIE is less singular than the kernel in 
the EFIE. As a consequence, less sophisticated expansion functions may be 
employed for representing the unknown current. On the other nand, the MFIE 
becomes unstable for thin structures. The vector cross product between J and R in 
the Green’s function may lead to computational difficulties by virtue of the small 
angle subtended, as in the case of the thin-wire [Ref. 9]. There is one unique feature 
of these integral equations for the induced surface current, J s , which is important to 
the discussion of natural resonance scattering. The unknown surface current, J s , 
inside the integral in both the MFIE and EFIE has the retarded time r argument. 
The retarded time r=t-R/c is always less then t since the P.V. integi .j ^ciudes the 
point R=0. Considering this feature, the unknown current J s (t), at any given time 
t, can be calculated from the MFIE in (2) using the known incident field at time t 
combined with the integral of terms which are known from the past history of the 
current. From another point of view, the effect of the source current at any point r' 
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is delayed by a time R/c in affecting the current at the observation point r. This 
point forms the basis of an iterative technique, termed time-stepping, for constructing 
the solution. The surface current may be determined by "stepping on" in time, thus 
eliminating the matrix inversion required for the numerical u.v.e of the frequency 
domain integral equation. Figure 2 shows the region o. • i* 'au? space-time for 
the one-dimensional case (x,t). This region is defined by ct-l .-x’ 1 <0 and denoted 
by the shaded area. 


Ct 



Figure 2. Region of Interaction 

The same general procedure can be similarly employed for solving the EFIE 
(1). The current at the spatial point of observation, anc at some time t, can be 
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calculated in terms of the known incident electric field at that same spatial point and 
time, and from the integral of terms involving retarded time [Ref. 8]. This point is 
further illustrated in the case of a thin-wire. 

C. NATURAL RESONANCE SCATTERING 

The following discussion is based upon the theory described by Morgan [Ref. 
4j. The MFIE (2) describes the induced current on the surface of the object in terms 
of the incident field or "physical optics" part, and in terms of the "feedback" current. 
Figure 3 shows the situation in transient scattering. An incident field with short pulse 
duration illuminates the scatterer. In the case of a radar it can usually be considered 
a plane wave. Once illumination of the object is completed, for t>t D H inc =0, and 
only the source-free currents remain on the object, generating the natural modes. 
These modes are of the form J n (r) exp(s„l), where the natural resonance frequencies 
s n =CT n +j&) n are functions of the scatterer geometry and composition. The SEM 
representation of the current takes the form of a summation of damped sinusoids 
which can be expressed using complex exponentials as 

V.<) = £ , I >t L . 

«=-•» 

Since J(r',t) is real, the complex exponents s m , which are poles in the complex 
frequency plane, and the weight coefficients A m come in conjugate pairs. For a 
practical incident field, having a finite bandwidth, only a limited number (say N) of 
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natural resonance modes will be "significantly" excited. The scattered field which is 
generated by the induced current is composed of two parts; an early-time driven 
response and a late-time natural mode response [Ref 4]. The early-time scattered 
field can be described as a linear combination of two terms. The first is an aspect- 
dependent physical optics term. 



Figure 3. Short Pulse Wave Illumination [After Ref. 4] 


The second term describes the scattered field due to the source-free current 
distribution. These currents are integrated over the time-varying portion of the 
surface to get the contribution from all previously illuminated points which are dot- 
shaded in Fig. 3. This term can be represented by the SEM expansion having the 
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same exponential resonance terms with time-varying coefficients (SEM class 2). The 
late-time scattered field is due to the source-free currents that remain after the 
incident field has completed its illumination of the scatterer. The late-time starts 
after a delay of T e =T+2D/c from the initial impact of the incident field on the 
object, where T is the practical incident field pulse width, and D is the length of the 
object in the direction normal to the wavefront of the incident field. At time T 0 the 
integral is calculated over a fixed surface area, and thus the SEM representation for 
the late-time is a summation of the same exponential resonance terms with constant 
coefficients (SEM class 1). 

In summary, the monostatic transient scattered signal waveform can be 
described in the following form [Ref. 4] 

>(<)=>^0 •[ w-tK‘-T .„>] *yfi> w-Tj+m w 

where U(t) is the Heaviside unit step function, y E (t) is the early-time response, y L (t) 
is the late-time response, and N(t) describes the measured noise and other signal 
pollutants. The late-time portion of the scatterer response to such an incident field, 
with its unique features, is considered in the process of this research. The discrete 
form of the late-time response is obtained in order to set up the required ARMA 
model. 
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D. THIN WIRE CASE 


Consider the case of a perfectly conducting cylindrical thin-wire scatterer. The 
geometry of the problem is shown in Fig. 4. The radius of the wire, a, is small 
compared to the wire length, L. The radius, a, is also small compared to the shortest 
significant incident field wavelength, X. 


z 



Figure 4. Geometry of a Cylindrical Thin Wire Scatterer 

In this case the thin-wire approximation is applicable; that is, the azimuthal surface 
curre-'t is negligible compared to the axially-directed component. Then the surface 
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current density, J s , can be written as 


* 2na 


( 5 ) 


where 2 is a tangential unit vector pointing in the axial direction. The current is a 
function of z only. The incident electric field E inc has contribution to the surface 
current only in the 2 direction. Under these conditions the integral equation for this 
case can be written, using the EFIE type (1), as 


4nE^(z,T:)= 

i R at 


Uf^f-Ll z(z W)dr'dz' 

* J R 3 Ldz' 


O 0 


1 r (z -z!) d 
e o C o 


R 2 * 


IXz'rfdz' 


( 6 ) 


where R=r-r' is the vector from the current element along the z axis pointing to the 
point of observation on the surface of the wire, DRl=[(z-z’) 2 +a 2 ] 1/2 , while r=t-R/c 
is the retarded time. The integral in (6) is not a P.V. type integral since r' *r. 
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There are several approaches to solve the thin-wire integral equation arising in 
time-domain scattering problems. The interpolation procedure and the finite- 
difference approach are presented in [Ref. 8]. The general interpolation procedure 
is to subsectionalize the thin-wire by dividing it into N segments and then to define 
a set of basis functions for expressing the unknown current, I, in each of these 
segments. Similar segmentation is also necessary for the time domain, choosing the 
appropriate time interval with regard of the spatial interval. An interpolation scheme 
in time and space is then used to express the current at one node (space and time) 
in terms of the current values in the neighboring nodes. The final step is to describe 
the thin-wire integral equation using the expansions and apply point matching to 
generate the desired matrix equation. 

A second approach is to use the finite-difference method to approximate the 
differential operators appearing in the time-domain integro-differential equation. 
This approach was introduced by Sayre and Harrington [Ref. 10] where the EFIE is 
written in terms of the vector potential A. The solution for the vector potential A can 
be based on the finite-difference method, as applied to the driven wave equation. 
In fact, aside from the specified boundary conditions, the solution has the same form 
as will be employed in the acoustic string case in the following chapter. 

£. SUMMARY 

The late-time portion of the scatterer response can be represented by a 
weighted summation of natural modes. The modes are functions of the scatterer 
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geometry and not of the incident excitation. The EFIE which describes the surface 
current on the scatterer can be solved numerically for the case of a thin-wire. 

Although, in principle, the ARMA representation can be derived by discretizing 
an integral equation, the full topological connectivity wrought by the Green's 
function yields an imposing difficulty. A better approach is found by employing the 
finite-difference approach, with its "nearest-neighbor" connectivity. This will be 
presented in Chapter III. 
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III. ELECTROMAGNETIC AND ACOUS 1 Z SCATTERING EXAMPLES 

A. INTRODUCTION 

Open region electromagnetic and acoustic scattering and radiation problems can 
be formulated using the integral equation approach. Other techniques may also be 
applied to these problems, both in the frequency- and time-domains. One of these 
methods is the finite-difference scheme which provides a convenient means for 
deriving a time-stepping algorithm for solving the EFIE. In this work, an exclusive 
use of this method has been utilized to formulate the problem of transient scattering 
via the time-space wave equation, and to demonstrate the new approach of 
presenting the late-time scattering response in a discrete ARMA model. 

The objective of this chapter is to describe the numerical solution chosen for 
this research. The derivation of the EFIE expressed in terms of the vector potential 
is presented. An analogous case of a vibrating string of fixed length is then described 
in conjunction w : th various methods of solution to complete the theory of the analog 
form of the problem. A discrete form of the time-space wave equation using the 
finite-difference method is then presented along with pertinent numerical 
considerations. 

The wave equation with forcing term describes numerous physical phenomena 
such a driven finite string or an illuminated transmission line. A computer program 
entitled TH7.FOR was written to support this research, and a source listing is given 
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in Appendix A. The program provides the solution of the inhomogeneous time-space 
wave equation. The string (or transmission line) can be excited by a Gaussian pulse 
from different angles of incidence. The amplitude and width of the Gaussian pulse 
excitation, as well as the number of segments on the string can be independently 
selected by the user. Results in time-space are presented for different cases. Fast 
Fourier transforms (FFT) are used to obtain the frequency-domain results from the 
time-domain data. Results in these cases are also presented showing agreement with 
the basis of the natural resonance scattering theory. Time-domain data provided by 
the computer program was also used to check the ARMA model and is presented in 
Chapter IV. 

B. DERIVATION OF THE VECTOR POTENTIAL EQUATION 

The Electrc Field Integral Equation (EFIE) may be written in terms of the 
vector potential A. The derivation begins with a slightly different form of the EFIE. 
The scattered electric field, produced by the induced current, may be written in terms 
of the scalar potential $, and the vector potential A, as 

«r,0=~'V4>. ( 7 ) 

ot 

The scalar potential </> can be expressed in terms of A, via the Lorentz condition, as 

V*A(r,f) + p o s <p(r,t) = 0 . (8) 

ot 
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Applying the boundary condition of the electric field on the surface of a perfectly 
conducting body in order to write Eq. (7) in terms of the scattered field, and using 
the Lorentz condition to write Eq. (7) in terms of the vector potential, the tangential 
component of Eq. (7) can be written as 

t 

nxE Ue =nx—~ — f nxVV-Adt' , reS . ( 9 ) 

dt c lJ 


Equation (9) is differentiated with respect to time to eliminate the integral. This 
yields the following equation 

„ dE ** c » c?A 1 » . p n\ 

nx - =nx -nxVV'A, reS, 

dt 3t 2 C 2 

where the expression for the vector potential A in terms of the induced surface 
current J s is 



R 


( 11 ) 


Equation (10) may be reduced to the one-dimensional case of the cylindrical thin wire 
along the x-direction, with length L and radius a, under the assumptions described 
in Chapter II. The following equation is obtained 


d 2 A 1 d 2 A . 3E 
-= -4rtc - 

0x 2 c 2 dt 2 ° dt 


( 12 ) 
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The vector potential in Eq. (12) is defined in a slightly different manner from Eq. 
(11). The vector potential A is defined as 


L 





(13) 


with R=[(x-x') 2 +a 2 ] ly2 . 

For the purpose of this work it was found convenient to modify equation (12). 
The term which includes the incident field is replaced by the function f(x,t) which 
describes the pulse excitation. A term which describes loss was added in order to 
investigate cases where damping was included as well. The final equation, which has 
been solved numerically by computer program, has the following form: 


d 2 U(x,t ) 1 &U(x,t) t dU(x,t) 

dx 2 c 2 dt 1 * dt /( % 


where the function U(x,t) satisfies the expression, and £ is the positive coefficient of 
the loss term. When £=0 the equation is reduced to the lossless case which is 
simply the nonhomogeneous wave equation. The last step is to define the boundary 
and the initial conditions in order to completely describe this problem. 
Homogeneous boundary conditions have been set in this case to present a total 
reflection at the ends of a string. This formulation also represents the cases of short 
or open circuits at the ends of an excited TEM mode transmission line. In the 
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electromagnetic case the potential on the ends of the single wire scatterer is non-zero 
and is thus not properly represented by the wave equation solution considered here. 

In terms of the function U(x,t) the boundary conditions may be written as 
U(0,t)=U(L,t)=0. The initial conditions can describe the voltage or current on the 
transmission line before the incident field impacts. There is assumed to be no initial 
voltage or current. Likewise, the string is initially at rest in the acoustic case. These 
conditions are formulated by homogeneous initial conditions as U(x,0)=0, and 
dU(x,0)/dt=0, where 0< x <L. Note that homogeneous initial conditions may be 
replaced by some functions U(x,0)=h(x) and 5U(x,0)/3t=g(x), where 0< x <L, in 
order to make the partial differential equation (PDE) homogeneous. This point is 
further illustrated by the ARMA model presentation. 

C. VIBRATING STRING WITH FORCING FUNCTION 

The case of a thin wire is analogous to the problem of a vibrating string of fixed 
length. This problem is thoroughly discussed in the literature, including [Ref. 10] 
which covers this topic in considerable detail. The formulation of a vibrating string 
problem is 

i #vm a5) 

at 2 c 2 a 2 

with the boundary conditions U(0,t)=0, U(L,t)=0, and initial conditions U(x,0)=h(x), 
and dU(x,0)/5t=g(x). Several methods are available to solve the problem of a 
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vibrating string of fixed length. One way to solve it is by separation of variables. 
This yields the eigenvalues and eigenfunctions corresponding to the homogeneous 
boundary conditions. A Fourier series method is then used to satisfy the initial 
conditions. This is found to be a convenient means to solve the problem in many 
cases. However, difficulties may arise due to the complexity of some initial condition 
functions which must be integrated. One applicable example in the case of transient 
scattering could be the function which describes a Gaussian pulse excitation. 

An equivalent technique, but in some ways more useful, is the "method of 
characteristics". In this approach the d’Alembert’s solution is applicable [Ref. 11], 
with the following form 

x+ct 

u(x t t )= o, ± j g(x , )dx ,. do 

X — f 


The solution is valid for the region defined by 0< x-ct <L and 0< x+ct <L . This 
region is shown dot-shaded in Fig. 5 which describes points of position and time such 
that signals from the boundary have not already arrived. The modification to the 
solution is made considering the boundary conditions which, in turn, imply multiple 
reflections as illustrated in Fig. 6. The simplest way to obtain the solution is to 
extend the initial condition functions as odd functions (around x=0) with period 2L. 
With these odd periodic initial conditions, the method of characteristics can be 
utilized as well as d’Alembert’s solution (16). 
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D. FINITE-DIFFERENCE APPROXIMATION TO THE WAVE EQUATION 


The discrete form of Eq. (14) is required for the numerical solution. In general, 
numerical solutions are required in cases where the shape of the area of integration 
or changes to the boundary and initial conditions make analytical solutions 
impossible. These changes do not fundamentally affect finite-difference methods 
although sometimes modifications to the methods are necessary. In this work the 
finite-difference method was used to approximate the derivatives in Eq. (14), and to 
make use of the time-stepping algorithm previously mentioned. 

The finite-difference method may be derived based on a Taylor’s series 
approach. The idea is to approximate the function U(x) at a point near x=x 0 , e.g., 
(x 0 ±Ax), using the polynomial approximation of the Taylor series. Through the use 
of Taylor series, it is possible to approximate derivatives in various ways. A finite- 
difference approximation for a derivative can be written using forward difference, 
backward difference or central difference. The central difference is found to be more 
accurate [Ref. 11] and it was used in this work. 

Using the central difference formula, the first partial derivative is written as 

3U,„ s U(X'+Hx,yJ-U(X'-bx,yJ n71 

dx (W 2tTx ’ 

and the second derivative as 


&U, N U(x+Ax,y)-2U(x,y)+U(x-Ax,y) 
_ (Vo) --- 


( 18 ) 
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Figure 7 shows the discrete space-time plane of the finite string problem. The 
string length L is subdivided into N segments with Ax=L/N being the spatial interval. 
The time interval is chosen to be At==Ax/c, where c is the velocity of propagation. 
The notation U/, with the spatial segment i and the time segment j, is used to 
describe the function U(x,t), where x=iAx and t=jAt. 




i ■ • Li^J 


0 1 23 ... 


L 

N 


x.= iAx 


Figure 7. Space-Time Discretization 


The approximated PDE for Eq. (14) may be written as 




Ax 7 


At 7 


-?■ 


2A t 


=f( - 


(19) 
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The relation between the space interval and the time interval may be used in 
Eq. (19) to yield Eq. (20). Equation (20) is written for U/ +1 in a form called the star 
operator. The star operator is illustrated in Fig. 8. 



' 

U(i.j+1) 

/ 


-•i.j) 

\ 

✓ / 

N 

JLJ(ill.j) 

/ 

\ ' 

\ 

f(i.j) 

✓ 



✓ 

U(i.j-l) 




Figure 8. Star Operator 

Fquation (20) describes the standard five-point finite-difference approximation to the 
Laplacian V 2 with slightly different weights. 

ur l =A[vix + vlJfPUl" +Dfl, ( 2 °) 

where 

p=-^L-A, D=-AAx 2 , cl=cf. Ax . 

2+cl 2+cJ 

Note that when £=0, which is the lossless case, cl=0, A=l, P=-l, and D=-Ax 2 . The 
star operator is used to describe the finite-difference equation in a simple way. Note 
that in this case the value of U, J does not contribute since At=Ax/c. Given the initial 
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and boundary conditions, which are homogeneous in this case, the values of the 
function U/ at any point (i,j) can be calculated by marching forward in time using the 
star operator. This is the solution procedure applied in the computer program given 
in Appendix A. 

E. NUMERICAL CONSIDERATIONS 

In any application of numerical solutions, the questions of accuracy, stability, 
and information bandwidth must be resolved. This section includes a discussion of 
some numerical aspects which are relevant in this case. 

The finite-difference approximation to the first derivative is consistent, meaning 
that the truncation error vanishes as Ax-*0. Hence, the exact solution is expected to 
converge as the number of segments increases. The sampling rates in the spatial and 
temporal domains are key factors in determining the accuracy of the numerical 
solution. The spatial sampling rate should be high enough to adequately resolve the 
spatial variation of the incident field as it propagates past the scattered The time 
sampling space should be high enough to adequately resolve the time variation of the 
pulse excitation. However, the sample points in time are not indeoendent of the 
space interval. Correlation between them is required because of equivalence between 
space and time in the retardation effect. Since the interactions between the currents 
on different points on the scatterer depend upon the velocity of propagation, the time 
sample spacing, At, must be related to the space sample interval, Ax by cAt< Ax [Ref. 
9]. This is also known as the Courant stability condition for the wave equation [Ref 
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11]. The inequality is equivalent to requiring that the space sample points be at least 
as far apart as the distance the electromagnetic wave travels, with velocity c, in the 
interval between two sample points in time. The relation between At and Ax 
determine the stability of the solution. The idea is that the numerical process should 
limit the amplification of the initial conditions. In this simple case of a finite string, 
the time sample spacing is related to the space sample interval by cAt=Ax. 

Among the factors which determine the sampling rate are the shape and the 
width of the incident field pulse, the scatterer size relative to the pulse width, and the 
highest frequency natural mode to be obtained by the solution. A delta function 
space-time impulse whose frequency spectrum extends from zero to infinity with 
uniform amplitude is desired. However, this is impossible from a practical 
standpoint. The approximation is made by a Gaussian impulse since it rapidly decays 
to zero. The same property is applicable in the frequency-domain, where the 
amplitude rapidly decreases with increasing frequency. In order to adequately resolve 
the incident field in time and space the appropriate sample spacing, which result in 
reasonably accurate numerical results, have been found to be on the order of one- 
fifth and up to one-tenth the pulse width in time and space, respectively. The pulse 
width is determined by the scatterer size and the highest frequency information 
required. 







F. COMPUTED RESULTS 


A computer program was written to obtain time-domain results to support this 
investigation. The program numerically solves Eq. (14) using the finite-difference 
method. The solution is determined using the procedure described in previous 
sections and the pulse shape is Gaussian. Gaussian pulse parameters such as pulse 
width, amplitude, and angle of incidence can be changed independently by the user 
to provide results for different cases. Additional parameters such as string length, 
number of segments on the string, time delay for peak pulse impact, and the number 
of time steps to be computed can also be varied. The coefficient for the case of 
loss, can be altered by the user within the program to investigate lossy cases. Note 
that for $=0 the program will solve for the lossless case. Figure 9 shows the 
Gaussian pulse excitation for different angles of incidence. Broadside pulse excitation 
is shown in Fig. 9(a), and Fig. 9(b) shows the pulse excitation for a 60 degree incident 
angle. In both cases there are 10 time samples in the pulse width, where the pulse 
width is defined between the points at which the amplitude is 10% of the maximum. 

Figures 10 and 11 show some results in the space-time domain. Figure 10 
shows the displacement on the string along 200 time steps for the case of broadside 
excitation. A one-meter "electromagnetic string" was assumed, with €=3x10* m/s. 
The one-meter length was subdivided into 15 segments which results in a space 
interval of 1/15 meter. The time interval is related to the space interval by the 
velocity of light ,c, which gives At=(Ax/c)= 0.22 ns. The pulse excitation has a 
Gaussian shape with pulse width of 2.0 ns. Figure 10a shows the result for the 
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(b) 60 Degrees Excitation 
Figure 9. Gaussian Pulse Excitation 
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lossless case and Fig. 10b is the case with loss added to the system by setting the 
value of the coefficient $ in equation (20) to be non-zero and positive. In these 
cases, when the string is excited by a broadside Gaussian pulse, the displacement is 
symmetric along the string. Symmetry is observed about the center of the string due 
to the symmetric boundary conditions, initial conditions and driver. The difference 
between the lossless case and the lossy case is in the amplitude. The displacement 
of each segment is a periodic function with a period of 2N time steps where N is the 
number of segments on the string, except in the early-time when the incident pulse 
is still exciting the string. This point is illustrated in Fig. 11. In this figure, the 
displacement of segments number 2 and 8 are presented for the lossy case having 
periodic propagation except when the string is excited by the incident pulse. The 
period is 30 time steps which is exactly 2N for N=15 segments. The early-time is 
about 15 time steps including several time steps before the pulse starts. 

Figures 12a and 12b show results for the case ol 30 degree incident angle on 
the Gaussian pulse. All other parameters are the same as for the case illustrated in 
Fig. 10. As expected the displacement of the string is asymmetric. The early-time 
in this case is longer than in the case of broadside excitation since it takes more time 
for the incident pulse to complete its excitation of the string. The time-domain 
characteristics in these results did not change. Both have the same period (2N time 
steps). Figure 12a shows the results for the lossless case and Fig. 12b for the lossy 
case. 
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(b) Lossy Case 

Figure 10. String Displacement for Broadside Excitation 
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(a) Lossless Case 



(b) Lossy Case 


Figure 12. String Displacement for 30 Degree Incident Angle 
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Time-domain results were further processed to obtain frequency-domain results 
using the Fast Fourier transform (FFT). The objective was to check the time-domain 
results against the natural resonance mode theory. The frequency spectrum of the 
displacement on each segment was checked in order to analyze the resonance mode 
presentation as a function of the position on the string. The procedure includes 
computation of the spectrum of the displacement as a function of time for each 
segment on the string using the FFT. A representative example is shown in Figures 
13 and 14, where time-domain data of Fig. 12 was used in this case. There are 15 
segments on the string and the incident angle is 30 degrees. Figure 13 shows the 
frequency spectrum of the displacement of segments number 3, 6, and 11. It is clear 
that the modes have the same frequencies for these segments, although each with 
different energy depending on position. Figure 14 shows similar results to that of Fig. 
13, but for all the segments on the string. At this point a tentative conclusion may 
be made that the resonance mode frequencies are independent of the position on the 
string. Note that the FFT of 256 points was taken using the data including the short 
early-time data. In this case, the early-time data may be used since this portion also 
includes information about the modes. Figure 15 shows the frequency spectrum in 
the case of broadside excitation and with the same parameters as in the case of 30 
degree incident angle (Fig. 13). The differences between the results are explicitly 
presented. In the case of broadside excitation, only even modes are excited while in 
the case of the 30 degree incident angle both even and odd modes are excited. 
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Relative Power 



Figure 13. Resonant Frequencies (30 Degrees) 
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Figure 15. Modes of All Segments (Broadside) 
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The fact that these modes are aspect-independent may be checked in the 
frequency-domain by using the time-domain data for different incident angies. 
Figures 16 and 17 show results for various incident angles. The string was subdivided 
into N=25 segments. The Gaussian pulse width was set to 2 nanoseconds in order 
to allow 15 samples in the 10% pulse width. These parameters result in a sampling 
frequency to be 7.5 Ghz, thus providing a frequency limit of 3.75 Ghz, as determined 
by the Nyquist sampling theorem. However, the incident pulse bandwidth is less than 
1 Ghz, hence the modes are presented only within this bandwidth. The string was 
excited from incident angles of 0, 25, 50, and 75 degrees. Figures 16 and 17 show 
results of a 256 point FFT of the time-domain data for the same segment in each 
figure. These figures show results obtained on segment number 22 and 6, 
respectively. In both cases the modes excited for different incident angles appear to 
be the same. Note that for the 0 degree incident angle (solid line), only the even 
modes are excited. 
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Figure 16. Modes on Segment 22 
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Figure 17. Modes on Segment 6 












IV. THE DISCRETE ARMA MODEL 


A. INTRODUCTION 

Based on the theory of natural resonance scattering, a radar target can be 
considered as a Linear Time-Invariant (LTI) system. There are several ways of 
describing LTI systems. A linear system can be described by its impulse response, 
by means of linear difference equations, a system diagram or by the system transfer 
function with poles and zeros. The system transfer function model is extensively used 
to present natural resonance modes in transient scattering research. As previously 
mentioned, these modes are functions of the scatterer geometry and composition, and 
are independent of the incident excitation. A finite subset of the scatterer poles can 
possibly be used to represent the scatterer in the process of discrimination. 

An alternative linear system model may be used to describe a scatterer in order 
to find its natural modes in a relatively simple way. This new approach describes a 
system by means of difference equations. Different scatterers, which are considered 
LTI systems, may be described by means of linear constant-coefficient difference 
equations. Once such a model has been set up for a specific scatterer, the natural 
modes can be directly determined by the coefficients of the differential equation. 

In this work such a model has been set up for the late-time response of a 
vibrating string with forcing function. An alternate physical problem, which fits the 
mathematical model, is an illuminated TEM mode transmission line with either open 
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or shorted ends. In the equivalent electromagnetic case ihe late-time starts after the 
incident field has completed its illumination of the target, and is completely described 
using a weighted expansion of complex natural modes. During this period the 
scatterer acts as a recursive system, generating the natural resonance modes. The 
output of this system depends only on previous values of the output. This type of 
system can be described by an auto-regressive moving-average (ARMA) model. The 
coefficients for the recursive terms are constant since the feedback mechanisms on 
the scatterer are constant and time-independent. Moreover, the ARMA model which 
describes the system has the same recursive coefficients for all spatial points, as 
expected. 

Two cases are considered: lossless and lossy . The analytical solution and the 
numerical solution for the natural modes in both cases are presented including 
comparisons with computed results. Development of the ARMA model for each case 
is then presented. Demonstrations of these developments are included through 
examples for both the lossless and lossy cases. 

A variety of methods are currently used to estimate ARMA parameters. These 
methods have been used to estimates poles of scatterers by applying them to given 
data obtained by measuring the backscattering signal. In this work only the basic 
version of the Prony’s method is applied to data obtained by the numerical solution 
described in Chapter III. The objective was to verify the development of the ARMA 
model. Results of the ARMA model are compared with results obtained using the 
Prony’s method for pole estimation [Ref. 13]. 
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B. THE PROBLEM 


The problems of transient electromagnetic scattering or the equivalent problem 
of acoustic scattering may be treated in the sampled signal case by considering the 
scatterer as a linear system whose input and output satisfy a linear constant- 
coefficient difference equation of the form 

>.n)=X>, <21) 

*=1 *.0 

For transient electromagnetic scattering, it is assumed that no surface current is on 
the scatterer before the scatterer is excited by the incident field, while the scatterer 
is initially at rest in the acoustic case. Therefore, the system may be considered 
casual, linear, and time-invariant. The a's and the b's in this case are real constants 
and the difference equation in Eq. (21) can be used to compute the output 
recursively [Ref. 6]. Considering the late-time case, the input x(n) is zero for n>n 0 , 
where the discrete time n G corresponds to T 0 in the analog case. This is the time at 
•vhich the incident field has completed its illumination of the scatterer. Hence, the 
difference equation for the late-time has the form 

N 

y(nh^2 a ky( n ~ k ^ n>n o • ^ 

*=i 

Equation (22) describes the system model for the late-time where the unknowns are 
the coefficients a v a 2 ,..., a N , and the number of delays N. Figure 18 shows the form 
of the system diagram, which serves as a graphical way of representing the same 
information contained in the difference equation (22). Equation (22) is referred to 
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as the homogeneous part of the general form of an Nth-order linear constant- 
coefficient difference described by Eq. (21) with a 0 = 1. 



y{n) 


Figure 18. Realization of the General Model [After Ref. 5] 

The homogeneous equation (22) has a family of solutions of the form 

yW'IXC. (23) 

m»i 

where z m are complex numbers. A unique solution requires a set of N auxiliary 
conditions since th°re are N undetermined coefficients. Substituting Eq. (23) into Eq. 
(22), the complex lumbers z m must be roots of the polynomial 

EV-o, < 24 > 

*=o 

assuming that all N roots of the polynomial in Eq. (24) are distinct. Based on the 
theory of natural resonance scattering, these roots are the poles in the z-transform 
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as the homogeneous part of the general form of an Nth-order linear constant- 
coefficient difference described by Eq. (21) with a 0 = 1. 



y(”) 


Figure 18. Realization of the General Model [After Ref. 5] 

The homogeneous equation (22) has a family of solutions of the form 

y(") = IXC> (23) 

m«l 

where z m are complex numbers. A unique solution requires a set of N auxiliary 
conditions since th°re are N undetermined coefficients. Substituting Eq. (23) into Eq. 
(22), the complex lumbers z ni must be roots of the polynomial 

f> t z-*=0, (24) 

*=o 

assuming that all N roots of the polynomial in Eq. (24) are distinct. Based on the 
theory of natural resonance scattering, these roots are the poles in the z-transform 
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x and t may be separated by substituting the solution (26) into the PDE (25) to 
obtain the eigenvalue problem 


fix 2 



(27) 


with the boundary conditions w n (0)=« n (L)=0. The general solution for the 
eigenvalue problem has the form 

u n (x)=c 1 e s ** lc +c 2 e~ s * l/c . (28) 

The boundary conditions are applied to obtain the solutions for s n and the associated 
modes. For the lossless case s n =jti> n , where u> n =mrc/L for n=±l, ±2, .... Natural 
resonance modes have the form 


U n (x,t) =sin(~x)« ; “* f . (29) 

The trivial solution is obtained for n=0. The solution for the initial value problem 
may be obtained by writing the final solution as a superposition of natural modes 

oo 

VM- E A.U.W) , n* 0. ( 30 ) 

/|S-W 

The constants A n which are usually complex, can be found by applying the initial 
conditions and representing them by the appropriate Fourier series, then matching 
term by term. 
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An alternative way to find the natural mode expansion for the initial value 
problem is by using the fact that these modes appear in conjugate pairs for real initial 


values. In this case the solution can be represented by 

•o 

U(x,t) = £ c n sin cos (w„t+ <p n ), ( 31 ) 

/! = 1 ^ 

where U. n =U n \ The real coefficients c n and the unknown phase $ n are determined 
by the initial conditions using Fourier series to represent the initial condition 
functions and matching term by term. 

For the lossy case the solution represented by Eq. (26) may be applied, 
« n (x)=sin(mrx/L), but a slightly different solution is obtained for the modes. This 
solution is substituted into the lossy wave equation which has the form 

(32) 

a* 2 c 2 a 2 8f ’ 


where $>0 is the coefficient for he loss term. The eigenvalue problem for the 
spatial variable x has a different form, yielding a complex set of solutions for s n . 
Substituting Eq. (26) into Eq. (32) to separate the variables x and t, the following 
equation is obtained for the variable x: 


&u n (x) 

dx 2 


c 


S n =0 n + ^n‘ 


(33) 
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The solution for s n is found as roots of the quadratic equation 

Sn 2 +Zc\ + (.™)=0. (34) 

D. NUMERICAL SOLUTION FOR RESONANCE MODES: LOSSLESS CASE 
The solution for the natural resonance modes involves the finite difference 
discrete form of the wave equation, using the star operator. Figure 19 shows the 
space-time discrete domain as applied to this development. The finite difference 
approximation results in a discrete equation of evolution which may be written for the 
lossless case in the following form 

U(k) =A‘U(k- 1) - U(k- 2), ( 3S ) 

where 



u(x v t k ) 

u(x 2 ,t k ) 


0 1 0 ••• 0 

10 10 0 

II 

1 

u{x M ,t k ) 

II 

o •- ••• ••• 0 

0 0 10 1 

0 - 0 1 0. 


The vector U{ k) is composed of unknown nodal values at the k-th time step and A 
is an M x M sparse matrix with ones along the two diagonals adjacent to the primary 
diagonal, regardless of size. The number of unknown nodal values is M which is 
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related to the number of segments by M=N'-1, where N' is the number of 
segments. The discrete form of the natural mode solution has the separable form 




Z e v * 

fl 


where U n is a constant in k, and s n =jo n in the lossless case. 


(36) 


I- 



Figure 19. Space-Time Discretization for Mode Solution 
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Both u n and Z n should have acceptable values to solve Eq. (35). The solution for the 
natural modes may be obtained by substituting the solution in Eq. (36) into Eq. (35) 
to yield the following eigenvalue problem 


<37) 

Equation (37) represents M linearly-independent equations for the M unknown 
eigenvalues X n . After solving Eq. (37) for the eigenvalues, the modes can be solved 
by the relation between Z n and k n , via Eq. (37), as 


Z n + Z'n'=K ■ 
n n n 


(38) 


Equation (38) solves the natural modes. The number of modes is 2M since for every 
value of A. n there are two solutions for the Z n which are conjugate values. This can 
be explained by the fact that in the finite difference equation (35) there are 2M 
degrees of freedom since there are M unknowns, and each unknown requires two 
previous values. From this, it is concluded that 2M coefficients are expected to 
appear in the ARMA model. In terms of number of segments, the number of 
coefficients in the ARMA model is given by N=2(N'-1), where N' is the number of 
segments. 

The initial condition solution may be obtained by a superposition of the modes 
as 


m-i xw. 


71=1 


(39) 
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where N is the number of modes. The A„ terms are excitation dependent amplitudes 
and w n (k) are excitation independent modes. There must be at least N=2M modes 
to allow the superposition (39) to be complete for any given initial conditions. 

E. NUMERICAL SOLUTION FOR RESONANCE MODES: LOSSY CASE 
The same solution procedure followed in the lossless case is applied to the lossy 
case. The finite difference equation may be written as 

U(k) = c x A'U(k- 1 ) - cJj(k-2 ), ( 40 ) 


where c 1 =2/(2+A,c£Ax), C2=-Cj+(c£Ax)/(2+c£Ax), and both A and C/(k) have the 
same form as for the lossless case. Note that for the lossless case (5=0) Cj = l and 
c 2 =-l. The same solution (36) is applied with s n =cr n +jo> n . By substituting the 
solution (36) into the PDE (40) and separating the variables, the following eigenvalue 
problem is obtained 

cAU=(Z+c-Z:')U=X n U n . ( 41 ) 

I /i x n 2 n ' n n n 


Equation (41) represents M linearly-independent equations for M unknown 
eigenvalues. Solving Eq. (41) for the eigenvalues, the modes can be obtained by 
using the relation between A. n and Z n in Eq. (41) as 




( 42 ) 
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The initial value problem may be solved by a superposition of the natural 
modes as in the lossless case, with N=2M modes to allow the superposition (39) to 
be completed for any given initial conditions. 


F. ARMA MODEL DEVELOPMENT 

The numerical solution shows that the modes are spatially-independent, i.e. all 
nodes have the same Z n . Hence the ARMA model may have the following form 


N=2M 

U(.k)= £ a m U(k-m), 


m-l 


( 43 ) 


where the a m ’s are the same for all spatial nodes, even those next to the boundaries, 
and U(k) is the superposition of all U n -Z n k , for ne[-M,M]. Note that the number of 
the coefficients is N=2M, where M is the number of spatial nodes along the string 
excluding those on the boundaries. The unknown coefficients can be determined by 
applying the z-transform to the recursion equation (43), which yields 

v£='L a .v£""- (44) 

m=l 


Equation (44) may be written for the coefficients in the form 


N 


Z N V - ' rjN-m 


( 45 ) 


1 
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or in the form 


ryN „ ~n- 1 , ~n- 2 „ y „ _ n 

a i^n ~ U 2^n ~"‘~ a N-\^n U N~®‘ 


( 46 ) 


Factoring the left hand side of Eq. (46) into first-order terms gives 


(Z-Z L ) (Z-Z_,) (Z-Z 2 ) (Z-Z_ 2 ) —(Z-Z W )(Z-Z. M ) 


( 47 ) 


where Z n are the poles and Z n =Z. n *. The last step is to use the results of the natural 
modes in Eq. (47) generating the polynomial and comparing term by term to Eq. (46) 
to find the values of the coefficients. 

The validation of the ARMA model can be obtained by showing that Eq. (43) 
can be derived directly from the equation of evolution (35) by working backwards. 
This point is further demonstrated using an example for the lossless case. 


G. VALIDATION EXAMPLE: LOSSLESS CASE 

Consider the case of N' =4 segments of an undamped string. The number of 
unknowns is M=N'-1=3. The vector of the unknown nodal values at the k-th time 
step has the form 


m= 


«(v*) 

u(x 2> t k ) 

u(x 3 ,t k ) 
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The equation of evolution in this case has the form 


with the matrix 


U(k)=A‘U(k-l)-U(k-2), 




0 1 0 
1 0 1 
0 1 0 


( 48 ) 


The number of discrete natural modes is N=2M=6. The eigenvalue problem with 
the form of Eq. (37) is obtained by substituting the separable form (36) into the 
equation of evolution (48). The eigenvalue problem has the following form after 
rearranging the equation 




( 49 ) 


Substituting A and solving for the eigenvalues gives three solutions: A n =0, and 
X n =±V2. The solution for the modes is obtained via 

(50) 

The solution has six modes: Z n =e ±j " /2 , e ±j " /4 , e ±j3,r/4 . Figure 20 shows the six poles 
in the z-plane. The poles are on the unit circle since this system is lossless. Since 
Z n =e je , where 0=a> n At, the modes are the same as in the exact case 
« n =n7r/4 A t=n-jrc/L. 
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The ARMA model has the form 


6 

m = Y,a m U(k-rn). 


m=l 


The z-transform of Eq. (51) is 


( 51 ) 


Z*=£ 




(V-m) 


m~ 1 


The coefficients a m can be found by factoring 


( 52 ) 


(Z-Z 1 )(Z-Z,*)(Z-Z 2 )(Z-Z 2 *)(Z-Z 3 )(Z-Z 3 *), ( 53 > 

with the known solutions for the Z n . The polynomial obtained by this procedure is 

Z 6 +Z 4 +Z 2 +l, ( 54 ) 

which gives the coefficients: fl,=0, a 2 =-l, a 3 =0. a 4 =- 1, fl 5 =0, and a 6 =- 1. The 
ARMA model for this case has the form 


U(k) = - U(k-2) - U(k- 4 ) - U(k-6). ( 55 > 

Verification of the ARMA model may be done by using the finite-difference 
equation of evolution (48) to find the coefficients in Eq. (55). The procedure is to 
use Eq. (48) to write f/(k-l), U(k-3). and JJ{ k-5) to finally obtain the form in (55). 
This is shown in Appendix B. 

In a case of a 5-segment string, there are 4 unknowns for each time step, and 
therefore 8 modes. The poles obtained in this case are Z n =e~ ) ~ !S , e" j2 " /5 . e** 3 ^, and 
e *j4r/s poles are shown in Fig. (21) on the unit circle of the z-plane. 
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H. VALIDATION EXAMPLE: LOSSY CASE 


Consider the case of N'=4 segments on a damped string. The number of 
unknowns is M=N'-1=3. The vector of the unknown nodal values at the k-th time 
step has the form of that in the example of the lossless case. The equation of 
evolution in this case has the form 


where 


U(k) = c x A-U{k- 1) - c 2 U(k-2), 


2 c$Ax 

c, =-, c, = — 2 -c,. 

2+c$Ax 2+cijAx 


(56) 


The same matrix >4 as in Eq. (48) is applicable in this case. The number of discrete 
natural modes is N=2M=6. The eigenvalue problem with the form (37) is obtained 
by substituting the oeparable form (36) into the equation of evolution (56). The 
eigenvalue problem has the same form as in the lossless case. Substituting A and 
solving for the eigenvalues gives three solutions: A. n =0, and X n =±c 1 V f 2. The solution 
for the modes is obtained via 


^n + c 2^n ~^n ^n =C lfe’ C \'J 2. 

The solution has six poles, 


(57) 


, */'■£ ±\l2c , ±t/2cf+4c, 

xj-m* - z »- — ■ — 
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Figure 22 shows these six poles in the z-plane. The poles are inside the unit circle 
and on a circle representing equal loss for all modes. In the electromagnetic case, 
as the frequency increases the higher frequency modes have more loss. 


Z-Plane 



Real z 

Figure 22. Modes of 4 Segment Damped String 
The ARMA model has the form 

6 

( 58 > 

m=l 

The same procedure used to find the coefficients for the lossless case is applied in 
the lossy case. The coefficients a m can be found by factoring the solutions for Z n . 
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The coefficients in this case are: ^=0, a 2 = 3c 2 +2cj 2 , a 3 =0, a 4 =-3c 2 2 +2c 2 c 1 2 , a 5 =0, 
and a 6 = c/. Note that the non-zero a’s are negative numbers with decreasing 
absolute values as their index increases. The ARMA model in this case has the form 

U(k) = (3c 2 +2cf) U(k-2) ~(3c2+2c 2 cf) U(k- 4) + c 2 3 U(*-6). ( 59 > 

Checking the coefficients by using Cj=1, and c 2 =-l yields the same coefficients 
obtained in the lossless case. 

I. COMPUTED RESULTS 

The objective of this part of the work was to obtain additional verification for 
the ARMA model. Tbe idea was to apply an algorithm which can estimate ARMA 
model coefficients to the time domain data generated by the computer program 
TH7.FOR. There are several methods available to estimate these parameters. 
Among them is the Prony’s method [Ref. 13]. The Prony’s method implemented via 
a computer program entitled TEST.M. A source listing is given in Appendix C. The 
program implements the following procedure for each segment on the wire. 

Using the time-domain results of a given segment, the following M x M matrix 
and vector are generated 



■«(M) 

u(t-2) ■ 

u(t-M) 


u(t) 

p= 

u(t- 2) 

u(t- 3) • 

■■ 

, v= 

«(f-1) 


u(t-M) 


- u(t-2M-l) 


u(t-M+l) 
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where M is the number of coefficients, and t is the time step at which the ARMA 
model is applied. The time-step t should be "late" enough so the time step t-2M-l 
is in the late-time portion of the time-domain solution, u(t) of the given segment. The 
vector a is composed of the coefficients a v a 2 , ...a M . The a m ’s coefficients are 
computed by 

a=P~ l -V, < 60) 

where P 1 is the inverse matrix of P. This method was applied to various cases. The 
results show that the coefficients are the same for all segments. The coefficients are 
also the same for any time-step t in the late-time, t-2M-l>n 0 , where n G is the initial 
time-step of the late-time portion of the system response. For all lossless cases the 
results show coefficients with the form aj=0, a 2 =-l, a 3 =0, a 4 =-l, ...a M =-l. In all 
lossy cases the results are a,=0 for odd i, and decreasing absolute values of negative 
numbers a„ for even i. In some cases the results were different from segment to 
segment. In those cases where the results produced segment-dependent coefficients 
the frequency response showed that not all of the modes where excited. A 
representative result is shown in Fig. 23. In this case, the damped string was 
subdivided into 11 segments, hence the number of modes is 20. The same 
coefficients were obtained for all segments, for various incident angles, and at each 
time-step in the late-time. 
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VI. CONCLUSIONS 


The solution of electromagnetic natural resonance modes by means of finite 
difference equations promises to be useful and effective. The ability to construct 
discrete ARMA type models of complex scatterers may be applicable to estimating 
their electromagnetic signature. In this thesis a simple ARMA model for the case of 
the one-dimensional wave equation was investigated in order to examine this new 
approach. Complex three-dimensional scattering structures can be treated using a 
space-time finite-difference approach which is an extension to that investigated 
herein. The ARMA model was constructed for the late-time response when the 
scatterer acts as a recursive LTI system. The process of constructing such a model 
included two main steps. First the scatterer should be presented in a discrete form. 
In this work the finite- difference approach was presented as well as the time-domain 
integral equation. The finite-difference approach is found to be convenient for the 
construction of the required model. The advantage of this method is in the locally- 
connected discrete form which can be effectively used in the next step. In the second 
step, an algorithm is constructed to explicitly present the natural modes. In this 
work, the algorithm is based upon the star operator obtained from the finite 
difference approach. It was shown that the star operator was used recursively to 
yield the ARMA model, with constant coefficients. 
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Results of this work have shown again, in the case of acoustic scattering, that 
natural resonance modes of a scatterer are spatial- and aspect-independent. The 
scatterer acts as a complete recursive system only in the late-time. Since the early- 
time is important because of practical considerations of signal to noise ratio, further 
work is required to extend the model to and to examine its complexity in the early- 
time. The computer program entitled TH7.FOR can be modified to help this 
investigation. 

The electromagnetic thin-wire scattering case was formulated via the vector 
potential. However, the potential on the ends of the wire is unknown. Assuming 
homogeneous boundary conditions reduces the problem to that of a vibrating finite 
string with zero displacement at the ends. Extension to the electromagnetic case may 
be made by attempting to discretize the current on the wire using the EFIE. In the 
case of currents, the boundary conditions are zero. 

A better approach is found by employing the finite-difference approach. In this 
method the spatial point current is given in terms of its "nearest-neighbor" currents 
yielding the applicable star operator. The ARMA model is constructed using the star 
operator. 


66 





APPENDIX A. SPACE-TIME WAVE EQUATION PROGRAM 


The program entitled TH7.FOR numerically solves the displacement of a finite 
string. The equation for the position f/(x,t) obeys the wave equation including a term 
of loss 

d 2 U(x,t) 1 #U(x,t) R dU(x,t) f(xf) 
dx 2 c 2 dt 2 ar 

The program uses the finite-difference method, with central difference, to solve 
the partial differential equation (PDE). The string with length L meters is subdivided 
into N segments. The boundar' conditions £/(x=0,t) and U(x=L,t) are set to zero. 
The excitation is a Gaussian pulse with the form 

/(*,0 = O4-ex;p-M'-W-W 

The pulse width is set between "10%" points (T1 in the program), where Theta is the 
incident angle. The time-delay to peak pulse impact is T2. The times Tl, T2 and 
Theta can be defined in the program, along with the amplitude of the Gaussian pulse, 
GA. Both space and time intervals are defined as S and T, respectively. 

c VARIABLES DECLARATION 
c 

INTEGER N.NNEXT.Th 

REAL L,S,GA,C1 ,R,A.B.P.D.UI ,U2,Tmax,T1 ,T2,A2,A1 ,PL(403) 

REAL PI,Tht,T3 
CHARACTERS TITLE 
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DIMENSION U(103,403), F(103,403) 
c 

write(*,*) 
write(*,*)’ 
write(*,*) 
writej*,*)’ 
wrrte(*,*) 
write (*,*)’ 

write(V)' 

write (*,*)’ 
write(*,*)’ 
write(*,*)’ 
write (*,*)’ 
write(*,*)’ 

write (*,*)’ 
write(V)’ 
write(*,*)’ 
write(*,*)’ 
write(V)’ 
write (*,*)’ 
write(*,*) 

write(**)' 

write(*,*) 
write(V) 
pause 
c 

c INITIALIZATION 
c 

WRITER,*) 

WRITE(*,*) ’Enter number of your MONITOR Type’ 
WRITE(V)' CGA = = = > Enter 0 ’ 

WRITE(V) ’ EGA ===> Enter 1 ’ 

READ (*,*) NS 
WRITER,*) 

WRITE(*,*) ’Enter number for 'Form Feed’ after plot’ 

WRITER,*)’ Enter 0 ===> NO FORMFEED (Laser)’ 

WRITE(V)’ Enter 1 = = = > FORMFEED (Impact)’ 

WRITE(*,*j 

READ (*,*) NFD 

ILINE=2 

Pi=3.14159 

C0=2.9979E+8 

Th'(=0.0 

Th=0 

a1=(-1.0)*LOG(0.1) 

WRITER,*) 

WRITE(V) ’ ENTER : String length in Meters’ 
WRITER,*) 

READ(*,*) L 
WRITE(V) 


PROGRAM TH7.FOR COHEN YUVAL MAY 1990 ’ 
SPACE-TIME WAVE EQUATION PROGRAM 

The program solves numerically the displacement’ 
of a finite string.The equation for U(x,t) is the ’ 
wave equation including a term of loss. 

The program uses the ‘Finite Difference Method’ ’ 
to solve the Partial Differential Equation (PDE) ’ 

The string with length L meter is suDdivided into * 

N segments. The boundary conditions are 
U(x,t)=U(L,t)=0 and initial conditions U(x,0)=0 ’ 
and dU(x,o)/dt=0. The excitation is a gaussian ’ 
pulse: G=Ga*EXP(-a*(t-Tmax-x*tan(Theta))**2) ’ 

The pulse width is set between ‘10%* points. The ’ 
incident angle is Theta. The space step is S and ’ 
the time step is T.’ 

Ready to begin ? 
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WRITE(*,*) 

WRITE(*,*) ' ENTER : Gaussian pulse parameters’ 
WRITE(V) ' Enter: Amplitude GA’ 

WRtTE(V) 

READ(*,*) GA 
WRITE(*,*) 

WRITE(V) 

WRITER/) 

Tl=1e-9 

T2=0.6*T1 

C 

CONTINUE 

NNEXT=0 

WRITER,*) ’ ENTER number of segments N’ 

WRITE(V) ’ (Any integer between 2 to 100)’ 

READ(V) N 
WRITE(V) 

S=L/N 

c * choose S/T=C* 

T=S/C0 

WRITE(V) ’ How many TIME steps to compute ?’ 
WRiTE(*,*) ’ Enter Integer between 2 to 400’ 

WRITE(V) 

READ(V) M 
WRITE(V) 

WRITE(V) ’ Enter: coefficient R ’ 

WRITE(V) ’ (R>0 or R=0) in IE-10 Units’ 

WRITE(V) ’ The case R=0 is LOSSLESS case.’ 
WRITE(V) 

READ(V) R 

R=R*1E-10 
Cl =C0*R*S 
A=2/(2+C1) 

P=(Cl/(2+Cl))-A 

D=A*(S**2)*(-1.0) 

IF(NNEXT.EQ.O) GO TO 8 
c 

5 CONTINUE 

WRITE(V) ’ Enter Pulse width in NSEC’ 

READ(V) T1 
Tl=Tl*1E-9 

WRITER,*) ’ Enter Time delay in NSEC for Peak Pulse Impact’ 
READ(V) T2 
T2=T2*1E-9 
NNEXT=0 
8 CONTINUE 
Tmax=T1/2 
a2=al/(Tmax**2) 

IF(T2.GT.Tmax) GO TO 9 

WRITE(*,*)' Time delay must be longer then half pulse width!' 
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TRY AGAIN !!!’ 


WRITE(V) ’ 

GOTO 5 

6 CONTINUE 

7 CONTINUE 

WRITE(* ,*) ’ Enter Incident Angle in DEGREE (0< =Theta=<90)’ 
READ(V) Th 

IF (Th .GT. 90 .or. Th .LT. 0) then 

WRITE(V)’ TRY AGAIN! 0=< Theta <=90 ’ 

GO TO 7 
END IF 
NNEXT=0 
9 CONTINUE 
c 

c GENERATE BOUNDARY and INITIAL CONDITIONS 
DO 300 j=1,M+2 
DO 400 i=1,N+1 
U(i,j)=0. 

400 CONTINUE 
300 CONTINUE 
c 

c GENERATE GAUSSIAN PULSE EXCITATION E=G*EXP(-a*(t-Tmax)**2) 
c PULSE WIDTH at '10%’ POINTS, 
c 

Tht=(Th/180.0)*Pi 
T3=T*tan(tht) 

DO 600 i=1,N 
DO 500 j=1,M+2 

F(i,]>GA*exp((-1.0)*A2*((j-1)*T-T2-((i-1)*T3))**2) 

500 CONTINUE 
600 CONTINUE 

WRITE(*,*) ’ Want to create file *f.mat* with driver data ?’ 

WRITE(V)' NO = = > Enter 0’ 

WRITE(V)' YES = = > Enter V 

READ(V) I 
IF (I .EQ. 0) GO TO 23 
OPEN(4,file=’f.mat’) 

DO 100 j=1,M+2 
DO 200 i=1,N 
WRITE(4,*) f(i.j) 

200 CONTINUE 
100 CONTINUE 
CLOSE (4) 

23 CONTINUE 
c 

c DISPLACEMENT COMPUTATION : MARCHING FORWARD IN TIME 
c 

DO 700 j=3,M+2 

jH*2 

DO 800 i=2,N 
ii=i-1 


70 



U1 =A*(U(i+1,j-1)+U(i-1,j-1)) 

U2=P*U(i,j-2)+D*F(ii,jj) 

U(i,j)=U1+U2 
800 CONTINUE 
700 CONTINUE 
c 

c CREATE a DATA FILE for OUTPUT POST PROCESSING (matlab) 
WRITER,*) ’ Want to create file ‘U.maf with displacement data ?’ 
WRITE(V)' NO ==> Enter 0‘ 

WRITE**,*)’ YES = = > Enter 1’ 

READ**,*) I 

IF (I .EQ. 0) GO TO 24 

OPEN(3,fiIe=’U.mat’) 

DO 710 j=2,M+2 
DO 720 i=1,N+1 
WRITE(3,*) U(i,j) 

720 CONTINUE 
710 CONTINUE 
CLOSE (3) 
c 

c OUTPUT 
c 

24 CONTINUE 

IF(NNEXT.EQ.O) GO TO 10 

2 CONTINUE 
WRITER,20) M 

20 FORMAT*’ Enter TIME in Time Step Units: [t=1,2...M=\l3,T) 
READ**,*) J 

OPEN(2,FILE='X') 

WRITE(2,11) J,N,Th 

11 FORMAT (’Displacement at TIME t= , I l3 r ’ STEPS. \i3, 

S’ SEG., Theta=’,i3) 

REWIND(2) 

READ (2,12) TITLE 

12 FORMAT(A) 

REWIND(2) 

NPTS=N+1 

Xmin=0 

Xmax=N 

j=j+2 

DO 900 i=1 , N+1 
PL(i)=U(i,j) 

900 CONTINUE 

CALL PLTSUB(title 1 npts,xmin,xmax 1 xmin I xmax t pl,ns 1 nfd 1 iline) 
GOTO 10 
c 

3 CONTINUE 
N1=N-1 

WRITE**,21) N1 

21 FORMATS Enter SEGMENT NUMBER: [1 I 2 I ...N-1=\I3 1 ’]’) 
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READ(V) i 
0PEN(2,FILE=’X') 

WRITE(2,13) i,N,Th 

13 FORMAT(’DISPLACEMENT of SEGMENT #*,13,’. \I3,* SEG., Theta=*,i3) 
REWIND(2) 

READ(2,14) TITLE 
REWIND(2) 

14 F0RMAT(A64) 
i=i+1 
NPTS=M 
Xmin=0 
Xmax=M 

DO 910 j=1,M+1 
PL(i)=U(i,j+1) 

910 CONTINUE 

CALL PLTSUB(title,npts,xmin,xmax,xmin.xmax.pl,ns,nfd.iline) 

GOTO 10 

4 CONTINUE 
N1=N-1 

WRITER,22) N1 

22 FORMATC Enter SEGMENT NUMBER: [1,2,...N-1 =*,13,*]*) 

READ(V) i 
OPEN(2,FILE=’X’) 

WRITE(2,15) i,T11,T22,N,Th 

15 FORMAT(’Driver on Seg.: , ,i3,*,Tpw= , ,F6.3, , ns,Td=*,F6.3, 

S' ns,N=’,l3, , I Theta=’,i3) 

REWIND(2) 

READ(2,16) TITLE 
REWIND(2) 

16 FORMAT(A64) 

NPTS=M 

Xmin=0 

Xmax=NPTS 

DO 920 j=1,NPTS 

PL(D = F(i*D 

920 CONTINUE 

CALL PLTSUB(title,npts,xmin,xmax,xmin,xmax.pl,ns,nfd,iline) 

10 CONTINUE 
WRITE(V) 

WRITE(V)' SELECT Number for Results or Change data* 


WRITE(V)’-’ 

WRITE(V) * Change Data AGAIN?? ===> 1* 

WRITE(V) * Displacement on the wire at time t=.. ===> 2* 

WRiTE(V) * Displacement on Segment #.. ===> 3* 

WRiTE(*,*) ’ Pulse Excitation as Function of Time ===> 4* 

WRITE( V)' Change Pulse Excitation Timing ===> 5* 

WRITE(V) * Change Pulse Excitation Incident ANGLE ===> 6* 

T11=T1*lE+9 
WRITE(V) 
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WRITE(*,17) Til 

17 FORMAT (' Current Values are: Pulse width=’,f6.3,’ ns.’) 
T22=T2*1E+9 

WRITER, 18) T22 

18 FORMAT (’ Time Delay for Peak Pulse impact=’,F6.3,’ ns.’) 
WRITE(M9)Th 

19 FORMAT (’ Incident Angle : Theta =\i3,’ Degree.’) 

WRITE(V) 

WRITER,*) ’ Any Ou .nr Integer = = = > END PROGRAM ! i!’ 
WRITE(*,*) 

READ(V) NNDCT 
IF(NNEXT.EQ.I) GO TO 1 
IF(NNEXT.EQ.2) GO TO 2 
IF(NMEXT.EQ.3) GO TO 3 
IF(NNEXT.EQ.4) GO TO 4 
IF(NNEXT.EQ.5) GO TO 5 
IF(NNEXT.EQ.6) GO TO 6 
STOP 
END 


c 

SUBROUTINE PLTSUB(nTLE,NPTS,XMil\!,XMAX,XS,XF,F,NS,NFD,ILINE) 
C 

C MS-FORTRAN Subroutine using "Grafmatic’ Library Subroutines. 

C Solid Line Using Portions of "PLOT* Program. 

C Written by M.A. Morgan with Latest Update to EGA/CGA June 1989. 

C 

C Default Printer is ’IBM Graphics* (e.g. Epson, Okidata, IBM) 

C With Plot Rotated 90 degrees From the Vertical. ■GrafPlus.Com’ 

C May be Run to Rotate Plot Upright on Paper and to Use a Variety 

C of Impact Printers. ’GrafLaser.Com* May be Run to Use a Laser 

C Printer. See GrafPlus/Laser Manual From Jewell Technology. 

C 

C INPU. DATA FORMAT 

C 

C TITLE - 64 Cha r acter Header 

C NPTS - # Data Points id 

C XMiN - Min X-Axis value 

C XMAX - Max X-Axis value 

C XS -Sturthg X value for F(1) 

C XF - Finisl ung X value for F(N) 

C F(N) - Input Data Array 

C NS - Monitor: *0’= GviA, *1*= EGA 

C NFD - *1* --> Form Feed After Plot Hardcopy (Impact) 

C Any Oth.rr Integer --> No Form Feed (Laser) 

C ILINE - *1 *--> + + + Symbol Dot Plot 

C Any Other Integer «> Solid Line Plot 

C 

C NOTE: The X-Axis Range Can Be Made Larger Than 

C The Domain of the Plotted Function, F(x). 
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Otherwise, Make XS=XMIN and XF=XMAX. 

CHARACTER*! DUM.BELL.FEED 
CHARACTER*64 TITLE.FNAME.HCOPY 
REAL X(512),F(512) 

INTEGER*2 N,JROW ) JCOL ) ISYM,ITYPE 1 NSCRN 
INTEGERS CYAN.WHITE.YELLOW.RED.BLACk.BLUE.NTWO 
INTEGER*2 JROW1,JROW2,JCOL1,JCOL2 
BELL=CHAR(7) 

FEED=CHAR(12) 

WRITER,*) BELL 
WHITE=7 
CYAN=11 
YELLOW=14 
RED=12 
BLACK=0 
BLUE=1 
NTWO=2 

NSCRN = 6 + 10*NS 
IF(NFD.EQ.I) ■ J P £N(1 .FILE^’PRN 1 ) 

N=NPTS 

DX=(XF-XS)/(N-1.0) 

FMIN=0.0 
FMAX=0.0 
DO 22 K=1,N 
X(K)=XS+(K-1.0)*DX 
IF(F(K).LT.FMIN) FMIN=F(K) 

22 IF(F(K).GT.FMAX) FMAX=F(K) 

IF(FMIN.GT.O.O) FMiN=0.0 
IF(FMAX.LT.O.O) FMAX=0.0 
C Computing Scale Factors for Vertical Axis 
ABSMIN=ABS(FMIN) 

ABSMAX=ABS(FMAX) 

YBIG=AMAX1 (ABSMIN.ABSMAX) 

N3CL=INT(LOG10(YBIG)) 

IF (YBIG.LT. 1.0) NSCL=NSCL-1 
YSCL=10.**NSCL 
FMIN=FMIN/YSCL 
FMAX=FMAX/YSCL 
ABSMIN=ABSMIN/YSCL 
ABSMAX=ABSMAX/YSCL 
DO 33 K=1,N 
33 F(K)=F(K)/YSCL 
C Adaptive Scaling Algorithm 5/89 
C Setting Polarity of YMIN and YMAX 
YMIN=0.0 
YBIG=0.0 

IF(FMIN.EQ.O.O) GO TO 37 
DY=0.5 
35 CONTINUE 
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IF(YMIh.GE.4.0) DY=1.0 
YMIN=YMIN+DY 
IF(ABSMIN.GT.YMIN) GO TO 35 
YBIG=YMIN 

YMIN=YMIN*FM1N/ABSMIN 
37 YMAX=0.0 

IF(FMAX.EQ.O.O) GO TO 41 
DY=0.5 

39 CONTINUE 

IF(YMAX.GE.4.0) DY=1.0 
YMAX=YMAX+DY 
!F(ABSMAX.GT.YMAX) GO TO 39 
IF(YMAX.GT.YMIN) YBIG=YMAX 
YMAX=YMAX*FMAX/ABSMAX 
41 CONTINUE 

C Calling GRAFMATIC Routines and Plotting F 
C Solid Line Graph Default 
ITYPE=1 
ISYM=-2 
NDOTS=0 

C + + + Line Graph If ILINE=1 
IF(ILINE.NE.I) GO TO 45 
ISYM=-1 
ITYPE=0 
45 CONTINUE 

JROW1 = 35+25*NS 
JROW2= 170+120*NS 
JCOL1 = 110-10*NS 
JCOL2= 540 
CALL QSMODE(NSCRN) 

CALL QPTXT(64,TITLE, YELLOW, 16,24) 

CALL QPLOT(JCOL1 ,JCOL2,JROW1 ,JROW2,XMIN,XMAX,YMIN,YMAX,XMIN, 
1 0.,0,1.,1.5) 

CALL QSETUP(NDOTS,CYAN,ISYM,CYAN) 

XMAJOR=(XMAX-XMIN)/5.0 

CALL QXAXIS(XMIN,XMAX,XMAJOR,1,1,2) 

IF(YBIG.LE.4.0) YMAJOR=0.5 
IF(YBIG.GE.5.0) YMAJOR=1.0 
IF(YBIG.EQ.8.0) YMAJOR=2.0 
IF(YBIG.EQ.9.0) YMAJOR=3.0 
IF(YBIG.EQ.10.) YMAJOR=2.0 
CALL QYAXIS(YMIN,YMAX,YMAJOR,1,1,1) 

JROW=32+21*NS+(ABS(YMIN)/(ABS(YMAX)+ABS(YMIN)))*(135+95*NS) 

JCOL=80-8*NS 

CALL QGTXT(3,’0.0',WHITE,JCOL,JROW,0) 

CALL QPTXT(1,’S’,YELLOW,5,18) 

CALL QPTXT(1 ,’c’,YELLOW,5,17) 

CALL QPTXT(1,’a’,YELLOW,5,16) 

CALL QPTXT(1,T,YELLOW,5,15) 

CALL QPTXT(1, ’e’,YELLOW,5,14) 
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CALL QPTXT{1, , X , ,YELLOW,5 1 12) 

CALL QPTXT(2,’10’,YELLOW,4,10) 

CALL QPTXT(1 YELLOW,5,9) 

CALL QPTXT(1 YELLOW,5,8) 

CALL QCMOV(0,8) 

WRITER, 150) NSCL 
CALL QTABL(rTYPE,N,X,F) 

C Hardcopy Query - Updated 5/11/89 
HCOPY=’Hardcopy ~> Enter P or p’ 

CALL QPTXT(30,HCOPY,RED,25,1) 

CALL QCMOV(55,1) 

READ(*,100) DUM 
HCOPY=’ 

CALL QPTXT(40,HCOPY,BLACK,25,1) 

IF(DUM.NE.’P’ .AND. DUM.NE.'p’) GO TO 48 
CALL QPSCRN 

IF(NFD.EQ.I) WRiTE(1,160) FEED 
48 CONTINUE 

C Exit to Blue Background on Screen - To Change This, 

C Replace ‘BLUE* in Calls to QPREG and QOVSCN to That Desired. 
CALL QSMODE(NTWO) 

CALL QPREG(O.BLUE) 

CALL QOVSCN(BLUE) 

100 FORMAT (A) 

120 F0RMAT(I5) 

130 FORMAT(E12.3) 

150 F0RMAT(4X,I3,\) 

160 FORMATC \A,\) 

RETURN 

END 
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APPENDIX B. ARMA MODEL ALGORITHM 


This Appendix presents the development of the ARMA model using the 
equation of evolution for the lossless case example given in Chapter IV. 

Start with the equation of evolution of the following form for M=3 

U(k) = A • U(k- 1) - U(k-2 ), W 

with the symmetric matrix 

0 1 O' 

A= 1 0 1 . 

0 1 0 

Using Eq. (1) to express U(k-l), substituting into Eq. (1) yields the following 

U(k) = - U(k- 2) +A 2 -U(k-2)-A-U(k~ 3). ( 2 ) 

Using Eq. (1) to express U(k-2) and substituting into Eq. (2) yields the following 

U(k) = - U(k- 2) +A 2 -[A-U(k~3) - U(k-4)] -A-Uik-3 ). ( 3 ) 

Rearranging Eq. (3) and using the fact that for M=3 A 3 =2A, the following 
expression is obtained: 

U(k) = - U(k-2) +A‘U(k-3) -£U{k-4 ). ( 4 ) 
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Using Eq. (1) to write U(k-3), substituting in Eq. (4), and adding and subtracting 
U(k-4) yields 


U(k) = - U(k- 2) -A‘U(k-5) + m-A) - U(k-4 ). ( 5 ) 


* 


Using Eq. (1) to write U(k-4) and substituting in Eq. (5) yields 


U(k) = - U(k-2) - U(k-4) - U(k-6). 


( 6 ) 


Equation (6) has the form of the ARMA model for M=3 described in Chapter IV 
Section G. 
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APPENDIX C. PRONGS METHOD PROGRAM 


The computer program entitled TEST.M implements the Prony’s method to 
estimate ARMA model coefficients using time-domain data. The program, which is 
written using MATLAB codes, finds the coefficients by implementing the procedure 
described in Chapter IV, Section I. The matrix with the time-domain data should be 
defined in a MATLAB format, and stored in a file called X.MAT. The data may be 
generated by TH7.FOR, and translated into MATLAB format using "translate" in 
MATLAB. The output of the program plots the results of the coefficients for each 
segment data. The coefficients are stored in a matrix called c where each column 
contains the coefficients for each segment data. Plots of results may be obtained by 
using the information in mat^'x c. 


% TEST.M program by Yuvai Cohen June 1990 
% The program implements the Prony’s method on time-domain data 
% estimating ARMA coefficients of the late-time solution of the 
% space-time wave equation. 

% 

load x 

input(’Enter number of segments of the string N‘) 
ns=ans; 

input('Enter number of coefficients M=2(N-1)’) 
m=ans; 

input(’Enter the time-step t’) 
t=ans; 

c=[l; 

for 1=1 :ns n=0; 

y=[J; 
u=[]; 
for i=1:m 
for j=l:m 
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y(i.i)=x(t-j-n,l); 

end 

n=n+1; 

end 

for i=1:m u(i)=x(t-i+1,l); 

end 

u=u’; 

a=y\u; 

c(:,l)=a; 

plot(a) 

end 
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